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ABSTRACT 

The interaction between settling of dust grains and magnetorotational instability 
(MRI) turbulence in protoplanetary disks is analyzed. We use a reduced system of 
coupled ordinary differential equations to represent the interaction between the diffu- 
sion of grains and the inhibition of the MRI. The coupled equations are styled on a 
Landau equation for the turbulence and a Fokker-Planck equation for the diffusion. 
The turbulence-grain interaction is probably most relevant near the outer edge of the 
disk's quiescent, or "dead" zone. Settling is most pronounced near the midplane, where 
a high dust concentration can self-consistently suppress the MRI. Under certain condi- 
tions, however, grains can reach high altitudes, a result of some observational interest. 
Finally, we show that the equilibrium solutions are linearly stable. 
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1 INTRODUCTION 

The presence of small dust grains in protostellar disks is 
critical to the thermal, dynamical, and chemical behavior 
of the gas, and is a crucial observational diagnostic. One 
particularly important feature is that dust tends to stabi- 
lize disks against the magnetorotational instability (MRI), 
and to otherwise complicate our understanding of MHD pro- 
cesses in such systems (ptone et al.||2000[ |Sano et al.||2000[ 



Salmeron fc Wardle|2008| |Bai fc Goodman|2009| ). By read 
ily adsorbing free electrons onto their surfaces the grains 
become charge carriers of very low mobility. In addition, the 
grains deplete the gas of alkali metal atoms, which are or- 
dinarily a low ionization potential source of electrons. This 
causes the resistivity of the gas-dust mixture to rise dramat- 
ically and essentially suppress MRI-powered turbulence over 
a wide range of heliocentric distances. Indeed, inclusion of 
dust grains in resistivity calculations result in greater "dead 
zones" than predicted by gas-phase chemical networks alone 
( Sano et al.|2000[|Bai fc Goodman|2009 i. 

Most studies thus far have assumed that dust grains 
were well-mixed over the vertical thickness of the disk (|Sano| 



et al. 2000 Salmeron & Wardle 2008 Bai & Goodman 2009 1 



This is justified for small grains, which are tightly coupled to 
the gas. However, as grains grow in size, while still control- 
ling the ionization of the gas, and gas becomes less dense, 
significant decoupling should occur, particularly at large he- 
liocentric distances. In particular, if the gas is stabilized, the 
embedded grains will settle toward the disk midplane. At 



this stage, depleted of its grains, the ionization rises and the 
gas may once again be vulnerable to the MRI. However, the 
ensuing turbulent agitation would resiir the grains, diffusing 
them once again upward into the gas (e.g. [Carballido et al!] 



|2006| ). Of course, the grains would then suppress the same 
instability that allowed them to diffuse through the gaseous 
envelope. The "cycle of inconsistency" would continue... 

How is this behavior ultimately resolved? It is possi- 
ble to envision a middle ground. In a turbulent medium, 
the growth of fluctuations is generally set by a balance be- 
tween linear (magnetic tension forces and resistive dissipa- 
tion) and nonlinear (cascade) processes. If saturation occurs 
at low amplitudes where nonlinearity is of secondary im- 
portance, the balance may be more simply regulated by a 
marginal, near zero, linear growth rate. The presence of some 
level of turbulent fluctuations would diffuse dust, and rather 
than cut-off the instability, such diffusion might regulate its 
growth. Specifically, the turbulence could stir just enough 
dust into the gas to ensure marginal growth: an increase 
in fluctuations raises the effective resistivity (more dust), 
a decrease in fluctuations increases the effective conductiv- 
ity (less dust). The question of the existence and stability 
of such a dynamical equilibrium is the focus of the current 
paper. 

Our chosen method of investigation is to construct and 
study a reduced system of ordinary differential equations, 
designed to reproduce certain key features of real disks. Our 
mathematical problem consists of two coupled, nonlinear 
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equations, both of which are a common staple of reduced 
systems. The first is a simple nonlinear Landau equation 
(Landau fc Lifshitz|1959 1 for the fluctuation amplitude. The 



second is a Fokker-Planck equation, with both drift and dif- 
fusive terms, for the concentration of dust grains. The diffu- 
sion coefficient of the Fokker-Planck equation is a function 
of the fluctuation level, whereas the growth rate of the Lan- 
dau equation depends on the grain concentration. It is this 
particular mathematical coupling in our proxy system that 
makes it interesting for astrophysical applications. We are 
able to demonstrate the existence of stable solutions for our 
reduced systems. 

This article is organized as follows: In Section 2, we in- 
vestigate the conditions under which the interaction studied 
is relevant in the disk. In Section 3, we outline our reduced 
system and its equilibrium solution. Proof of its stability is 
deferred to Appendix [B] In Section 4, we adopt a specific 
form for the growth rate for illustrative purposes. In Section 
5, we conclude. 



2 DUST AND MRI: A REVIEW 

In this section, we review the role of dust on the MRI in a 
schematic way to orient the reader with respect to order-of- 
magnitude scalings, and to highlight the conditions under 
which the interaction of dust settling with the MRI may be 
relevant. There are three primary criteria for this: 

(i) The dust dominates recombination of ions and elec- 
trons. 

(ii) Nonideal MHD effects are important, but do not sup- 
press MHD turbulence altogether. 

(iii) A substantial fraction of the dust grains can settle to 
the midplane. 

We shall quantify each of these in the next subsections. 

The disk is described in a cylindrical coordinate system, 
with R the heliocentric distance and z the height above the 
midplane. Since our calculations are local, we need not spec- 
ify a global disk model but we shall normalize our results 
to values of order those of the minimum mass solar nebula 
(MMSN; Hayashi||1981 \ near an heliocentric distance of 10 
AU. We will assume the disk to be axisymmetric and verti- 
cally isothermal, with the gas density given by: 






(1) 



with E the surface density, Cs the isothermal sound speed, 
H = Cs/ii is the pressure scale height and Q the Keplerian 
angular velocity. P = pc^ is the corresponding pressure. 



2.1 The ionization fraction 

Consider a gas composed of neutrals, ions and electrons, and 
a population of dust grains, of respective number densities 
Tin, Hi, Tie and Tip, all assumed to be at the same temperature 
T. The neutrals are predominantly H2 molecules and the 
ions are treated as one singly-charged species. The grains 
are assumed to be identical spheres of radius a and internal 
density ps ; we also denote by Pp = 47rpsa^Wp/3 the dust 



mass density and by e = pp/ p the dust-to-gas mass raticPj 
We denote by C, the ionization rate, which in the outer solar 



system shall be dominated by cosmic rays ( Bai fc Goodman] 
[2009| ). The evolution equations for Ue and rii are: 

Ik 



■ (rin ~ Pcsn^rii — I^-Ra vrenpUe 



(2) 
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where /3cff is the effective gas-phase electron-ion recombina- 
tion rate, VTe,i = \/8kBT/-jTmi^e 
speed, and le 



is a characteristic thermal 
the (averaged) product of the sticking coeffi- 
cient aB^__thejocusing^£tor due to electrostatic effects (the 
J of 



Drainefc Sutinll987l 






In the absence of dust grains, using charge neutrality, 
the fractional abundance of electrons at equilibrium ioniza- 
tion is ( |Gammie|1996l [Fromang et al.|2002| : 



Tie 



C 



/3cff"-n 



= 2x10-^° 



- C 10-^^mVslO-«kg/m^ \;:,^ 

^ 10-17 s-i /3eff p )^ ' 



Since metal ions have a low recombination rate (see Ap- 
pendix[A| , a small fraction of their cosmic abundance is then 
sufficient for Xe to allow widespread MRI activity ( Fromang 
|et al.||2002| |Bai fc Goodman||2009[ ) if dust is ignored. The 
appreciable depletion (>10 %) of chondrites in moderately 
volatile elements (e.g. alkalis) relative to the total condens- 
able matter (e.g. Scott fc Krot|2003 \ suggests that these ele- 
ments were not very efficiently removed from the gas phase. 
Their depletion in the gas phase likely did not exceed 1-2 
orders of magnitude then (similarly to cold interstellar gas, 
e.g. Yin||2005 I, in contrast to the much stronger depletions 



envisioned in some parameter studies (e.g. [Fromang et al 



2002[|Ilgner fc Nelson|2008[ [Flaig et al.|2012l ). Metals per se 
would thus be sufficiently abundant to significantly reduce 



the extent of the dead zone (Fromang et al. 2002 Bai fc 



Goodman|2009 1 but this does not hold if dust is taken into 



account ( Sano et al.|2000 Bai fc Goodman|2009 \ . Dust thus 
appears as the main agent acting to suppress the MRI. 

If gas-phase recombination can be neglected (which we 
seek here to quantify), we have, at equilibrium: 
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(5) 

The normalizing value of psa = 10^^ kg/m"^ applies to 
micron-sized grains. The smaller the grains, the larger the 
area offered for recombination per unit mass, and hence the 
lower the ionization fraction. Since the equilibrium attain- 
ment timescale a;e/C is shorter than all other timescales of 
interest, chemical equilibrium will be assumed throughout. 



1 The grains are considered to be sufficiently large to ignore the 
effects of their electrical charge on their own dynamics and their 
direct contribution to the current density ( see |WardIe|2007l l and 
thence nonideal MHD terms (but see |Bai| | |2011[ | for the effects 
PAH-sized grains). 
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Neglect of electron-ion recombination is warranted if: 
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(6) 



(Some justification of the normalizing value for /3eff is pro- 
vided in Appendix A, where the molecular ion/metal ion ra- 
tio is estimated). Under this condition, the dust may be said 
to control the ionization fraction as required by criterion (i). 
Note that our focus on dust properties is strictly justified if 
its impact on the ionization fraction is not overshadowed by 
vertical variations of the ionization rate, which is satisfied 
if e.g. the gas column density is smaller than the stopping 
grammage of the ionizing radiations (9.6 x 10^ kg/m^ for 
cosmic rays according tO|Umebayashi fc Nakano|1981 1. 



2.2 MRI activity 

The reduced ionization due to dust enhances nonideal terms 
in the induction equation, whose importance we now quan- 
tify. Various dimensionless numbers have been defined in 
the literature depending on the diffusivity regime, and the 
thresholds for MRI activation are still being debated, so we 
shall restrict ourselves to two of them, pertaining to ohmic 
and ambipolar diffusion, respectively: 

Ohmic diffusion is believed to dominate near the mid- 



plane (e.g. Salmeron & Wardle 20081. The importance of 



Ohmic diffusivity rjo is measured by the magnetic Reynolds 
number, which, if we inject equation (|5|, is given bjP] 
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m the neutral-electron cross-section 
(Draine et al. 1983 1. Currently estimated thresholds for Rei\/ 



for go od coupling between the gas and magnetic fields are 
10 



2±2 



(Fromang et al. 20021 



The importance of ambipolar diffusion, which may dom- 



inate in the upper layers of the disk ( Perez-Becker & Chiang 



2011 1, is measured by the dimensionless ion- neutral collision 



rate (per neutral molecule): 
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2 We use TjQ = AvXeTTLf^ 

and (25) of Balbusj (2011 1 



ap/i^'ipoe^x^) from equations (9), (24) 



with /3in = 1.9 X 10 ^^ m^/s the ion-neutral collision rate 
coefficient ( Draine et al.|1983 ) . The threshold for go od cou 



pling between ions and neutrals is of order 1 — 10 ( Perez- 
[i3ecker fc Chiang|2011||Bai fc Stone|2011| ) 



As mentioned previously, other dimensionless numbers 
exist in the literature but as these are proportional to the 
ionization fraction, their dependence on dust properties is 
subsumed in the factor psa/ pp, or equivalently psa/e. 

The condition (ii) that the turbulence level is a strong 
function of ionization is met when the magnetic coupling 
is marginally good, i.e. when the dimensionless number 
pertaining to the relevant diffusivity regime is close to its 
threshold value for ideal MHD. In other words, the region 
most relevant to our analysis is near the outer edge of the 
dead zone. From equations ([7| and (Isl, this edge may be 
expected to be at ~10 AU from the Sun in a MMSN, as it 



was in the detailed calculations of Bai & Goodman ( 2009 1 



2.3 Grain dynamics 

The dynamics of solid grains are dictated by gas drag. For 
grains small compared to the gas mean free path, the stop- 
ping time is (Epstein|1924l 



TT Psa 
8 pCs ' 



(9) 



For small grains [t <^Q, ^ ) , the vertical drift velocity of the 
grains is given by ( DubruUe et al.|1995 l: 

(10) 



v{z) = — f7 TZ. 

Settling is counteracted by turbulent diffusion, which tends 
to stir dust particles above the midplane. The vertical dif- 
fusion coefficient may be parameterized as: 



D, 



Q' 



(11) 



where 5z is a dimensionless parameter of order the standard 
a parameter (e.g. [Johansen et al.|2006[ ). 

With these notations, equilibrium between diffusion and 
settling is attained on a timescale 



f-vcrt — 
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(12) 



and the equ ilibrium thickness of the dust layer is Hp 
H/\/lT~S~z ( |Cuzzi et al.|l996[ ). We have introduced 

VLr 



_ TT psa 

^ 2 E(5^ 
= 20 



(13) 
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10-" 



102 kg/m2 



where r and 5^ are evaluated at the midplane {z = 0). S^ 
is thus a measure of the settling of dust relative to the gas. 
Significant settling of dust, with interesting feedback on the 
MRI, corresponds to Sz S> 1. This implies (1) relatively 
low surface densities (as expected far from the Sun) (2) rel- 
atively big grains, say 10^ pra (but likely no larger for 
dust to retain control of the ionization fraction), presum- 
ably as a result of coagulation or (3) a low turbulence level 
around the midplane, or a combination thereof. Low values 
oi 5z < \Q~^ are seen around the midplane in numerical 
simulations of layered accretion (e.g. Fleming fc Stone|2003 
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Ilgner fc Nelson|2008| |Oishi fc Mac Low|2009[ [Turner et al.| 



2010[|Okuzumi fc Hirose|2011||Flaig et ai.|2012[ ) . Note that 



while Turner et al. (20101 did account self-consistently for 



the dynamics of dust as well as its role on the ionization frac- 
tion, no significant effect of dust motion was found (see e.g. 
their figure 14), but this is because maximum (midplane) 
values of Sz were only 0.2, 0.8 and 1 for the runs with 1 /xm, 
10/im and 100 /xm grains, respectiveljrj 

With these numerical estimates setting the scales, we 
turn to a reduced model of interaction between dust and 
MRI turbulence. 



3.2 Equilibrium solution 

Equilibrium implies: 



d\ne v{z) 
dz Qy'^ 



y = 



A 



(17) 



(18) 



Plugging equation \1S\ in equation (llTl and integrating with 
respect to z yields: 



3 REDUCED MODEL 

3.1 General Equations 

Consider a vertical section of a protostellar disk. We denote 
by 1/ a turbulent fluctuation of the gas (with e still denoting 
the dust-to-gas ratio). Our reduced model consists of the 
following system of equations: 



de _ d 
^dt^ dz 



p[Qy 



dz 



':{z), 



g = ^{e)y - V 



Equation ( 15 1 may also be rewritten as: 



dt 



27(e)y" - 2Ay^ 



(14) 



(15) 



(16) 



The first equation ( 14 1 is a diffusion equation for the grains 
with a drift term proportional to the velocity v(z). It has a 



standard Fokker-Planck form. The second equation ( 15 1 is a 
Landau equation ( [Landau fc Lifshitz||1959"l ) for nonlinearly 
damped fluctuations. The constant Q relates y to the verti- 
cal diffusion coefficient D^ = Qy^ and the constant A char- 
acterizes the nonlinear saturation. Linear growth and damp- 
ing (here MRLdriven) are embodied in 7 the nominal rate 
coefficient, which depends on e by assumptions (i) and (ii) in 
Section 2. The form of 7(e) is at this point unprescribed, but 
it is assumed to be a positive, differentiable quantity. Phys- 
ically, we would expect it to be a monotonically decreasing 
function of e. It would reach an asymptotic value for e below 
the threshold for recombination on dust grains to be impor- 
tant (see equation (Im), which corresponds to the dust-free 
(not necessarily ideal) MHD turbulence value. For e above 
another threshold, MRI-powered turbulence is suppressed, 



but other insta bilities (e.g. (W cidcnschilling 1980; Youdin & 
|Goodman|2005| [Latter et al.|2010,|Lcsur fc Papaloizou|2010| ) 
may help maintain a minimum level of turbulence and es- 
tablish a dynamical equilibriumrl 



^ We identify 5^ with a evaluated from their figure 7. 

■* The Landau equation used here ignores any transport within 



the gas, e.g. from the active layers to the dead zone ( |Fleming| 
|fc Stone|[2003[ [Turner k. Sano||2008 ) but were we to add a term 



Jz (Qy'^-J^) in the Landau equation ||16| 

Ay"^ would be of order Qk'^/A ~ {n/^){S, 

k ~ Hp the reciprocal lengthscale of variation and it would 

thus be negligible. 



its ratio with e.g. 
+ l)&z < 1 with 



7(e) de 

(0) Idf e 



z-'/2H-' 



(19) 



where we have introduced 7^/ = 7(e = 0), the dust-free 
value of the growth rate, and Szdf = AT{0)c^/'ydf the corre- 
sponding value of Sz- The value of e(0) must satisfy: 



Pp(z)dz 



1 



2ttH 



e{z)e 



-z^ I2H' 



dz, (20) 



where we have introduced the dust-to-gas column density 
ratio e. If we approximate e(0) ~ e\/l + Sz (see Section 2.3), 
where Sz is evaluated at the midplane, this may be replaced 
by the following simplified relation: 



S: 



^d/— Z 



7d/ 



7(e^rTs;) ■ 



(21) 



In principle, depending on the mathematical expression of 
7(e) and on the values of Szoj and e, it is conceivable that 
more than one solution exists in terms of Sz , or, equivalently, 
e(0), or even none if 7 vanishes too rapidly with increasing 
e (which however seems unrealistic, see Section 3.1). 

We show in Appendix [B] that regardless of the form 
of 7(e), the equilibrium is linearly stable. From equation 



( B6 1 in the appendix, one estimate the damping timescale to 
be ivcrt(A///p)^ with A the vertical lengthscale of variation 
of the perturbation: this is basically a diffusion timescale. 
Physically, if we schematically distinguish between a "mid- 
plane zone" and an "atmosphere zone" , we may interpret the 
lack of instability due to a dust-controlled diffusivity as fol- 
lows: if the "atmosphere" has excess dust, diminishing tur- 
bulence there will make the dust fiow toward the midplane 
to cancel the corresponding dust depletion. If, on the other 
hand, the atmosphere has a dust depletion, enhanced turbu- 
lence will soak the excess dust from the midplane through 
the "interface" between the two. 



4 EXAMPLE 

As an example, consider the following functional dependence 
for the growth rate: 



7(e) = 



ydf 



1 + e/e. ■ 



(22) 



One may think of e* as the critical value of e for which 
the relevant dimensionless number of Section 2.2 is at its 
threshold value. The asymptotic value of the growth rate 
for e 2> e, is zero, i.e. we ignore any "background" hydro- 
dynamical turbulence. 
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Turbulence in a MM5N 




without settling 



self-consistent 



10^ 



Figure 1. Plot of the equilibrium dust-to-gas ratio e profile, as- 
suming the dependence of the growth rate on t in equation | |22| . 
Curves are drawn for S^^f = 1 (i.e., marginal settling for dust- 
free turbulence levels) for different values of e(0)/e* as marked on 
the figure, assuming a dust/gas column density ratio of 10^^. 



simplified self-consistent determination of settling 



Sj,Tar/T(e{I)» 




Figure 2. Plot of the left-hand-side (continuous) and the right- 
hand-side (dashed) of the simplified equation | |21| l for the example 
functional dependence of the growth rate displayed in equation 
([22b. For any parameter, there is always one unique solution (cor- 
responding to the intersection between the two curves) for this 
functional dependence. Here, e/e* = 10 and S^df = 1- 



We obtain from equation ( |19[ ): 
exp LSzd/(e 



1 1 



1) 



(23) 

This is plotted in Fig. IT] It is seen that enhanced dust abun- 
dance at the midplane entails an enhanced settling efficiency, 
hence a peaked distribution. But as 1^1 increases and thus the 
dust fraction drops, turbulence can increase and the profile 
becomes shallower, until dust no longer affects the turbu- 
lence level, by which altitude one has: 



,{z) 



exp( 



1 I 



e(0) 



1)) 



(24) 



that is, (1 -I- e(0)/e*) ^ times the result for a vertically con- 



Figure 3. Radial variation of the turbulence parameter 5z (eval- 
uated at the midplane) in a MMSN assuming that the growth 
rate is given by equation | |22| l. We assume a constant dust-to-gas 
column density ratio e = 0.01, paa = 0.1 kg/m^ (i.e. 30 /im- 
radius grains), S^af = 10~^ and take e, = 10~^(i?/10 AU)^/^ 
(to mimick the dependence on Q of dimensionless numbers per- 
taining to nonideal MHD in Section 2.2). The dotted line is the 
estimate ignoring vertical settling and the continuous line is the 
self-consistent estimate using simplified equation (|21|. 



stant 6z (e.g. Fromang fc Nelsonl|2009 |, a depletion due to 
the aforementioned enhanced dust concentration at the mid- 
plane. If e(0) ^ £,, the self-consistently determined dust 
concentration profile sets the vertical extent of the dead zone 
and the active layers. 

The mass balance constraint on dust expressed by equa- 



tion (20 1 reads here: 



2 



1 

<0) 



exp(5',d/(e-'-l))- 



-dx 



(25) 

Since the right-hand-side is a monotonic function of e(0), in- 
creasing from to +<X), there is always one unique solution, 
given S^df and e. 



The simplified equation (21 1, illustrated graphically in 



Fig. I2J may be solved for Sz as 

Szdf 



■>zdf 



1 + 



+ \ (l + Szdf) 



- \ 2 Q2 

ej 4 



- \ 4 



e» 



(26) 



If e ^ e, and Szdf(j/^*Y ^ li there is essentially 
no effect of dust and Sz ~ Szdf- At the other extreme, if 
S'zd/max [{e/e,), (e/e*)^) » 1, we have: 



■>zdf 



zdf- 



2 E5 



zdfC 



(27) 



where we have set e, = c^paa, with c« independent of dust 
properties, in order to account for the dependence of the 
ionization fraction on grain size (see Section 2.1). We see 
that, in this limit, the dependences of settling on grain size 
cancel out. This is because for a given dust concentration, 
larger grains drift more rapidly toward the midplane, but at 
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the same time allow higher turbulence levels. Certainly, this 
(asymptotically) exact cancellation is specific to the depen- 
dence we have chosen, but the simulations of [Turner et al.| 
(20101 show only a weak dependence of Sz on size, with 
the former varying by a factor of ~5 despite a two-order-of- 
magnitude variation of the latter (see Section 2.3). 

Extrapolating for a range of heliocentric distances, we 
have plotted the radial profile of S^ evaluated at the mid- 
plane in figure [3] It is seen that ignoring settling of grains 
leads to an overestimate of 5^ and hence an underestimate 
of the heliocentric distance of the outer edge of the dead 
zone. This is because settling induces larger dust concen- 
trations at the miplane, compared to the perfect vertical 
mixing assumption, and thus lower turbulence levels there. 
As to the vertically averaged S^ — a proxy for the standard 
a parameter — ignoring settling would lead to an underesti- 
mate of its value, since 7 is here a convex function of e, but 
in general, the effect of settling on this average depends on 
the particular mathematical form of 7(e). 



5 CONCLUSION 

We have studied the interaction between MRI turbulence 
and dust grains, allowing the dust to control the ionization 
level (and thus MHD turbulence level), while turbulent fiuc- 
tuations inhibit dust grain settling to the disk midplane. We 
have used a reduced model consisting of two coupled equa- 
tions: a Landau equation for the turbulent fluctuation am- 
plitude and a Fokker-Planck equation governing the vertical 
dynamics of the grains via a fluctuation-dependent diffusion 
coefficient. Unconditionally stable equilibrium solutions for 
the vertical grain distribution were found. Compared with 
models in which the turbulent fiuctuation level is constant 
with height, the solutions were signficantly more concen- 
trated near the midplane, with a lower level of MHD turbu- 
lence in this region. 

From simple estimates of the ionization fraction and 
the settling parameter, we found that the interaction stud- 
ied here is most relevant near the outer edge of the disk's 
dead zone. Indeed, the grains are likely to determine the 
dead zone morphology, and will likely control the region's 
heliocentric extent. Also, an enhanced settling could lead 
to flatter disks and hence steeper decreases of temperature 
with heliocentric distance. Thus, the model may be used in 
conjunction with numerical simulations (global or local), to 
sharpen observational predictions that will be of great in- 
terest when ALMA becomes fully operational. 
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APPENDIX A: ESTIMATION OF THE 
MOLECULE/METAL ION RATIO 

In this appendix, we provide some justification of the nor- 
malization chosen for P^a in Section 2.1. To that end, we es- 
timate the ratio between the number density n„+ of molecu- 
lar ions {m'^ , e.g. HCO'*') and that njv/+ of metal ions (M"*", 
e.g. Mg^, Na"""). We thus need to refine the model in the 
main text by distinguishing between the two. The resulting 
governing equations are the same as |Froniang et ar] ( |2002[ ), 
whose notation we adopt, with the addition of grains: 



chosen for j3cb in inequality m[. The latter then becomes: 

PcffniTie 



Cn„ 



^merrii Ij^j+ /3 i'tm+ P^a C 



mH2 hh Pt Xm e P 



10-^ 
lO'^kg/m^y \^ Xm 



10-17 s-iy V P 
Psa \ /lO"^' 



< 1 



(A9) 



™ : (n„ ~ I3nen^+ — j3tnMn^+ — /„,+ tto Vj,^+ npn^+ 



dt 



(Al) 



OTIm+ _ 7-2 

— -— = PtnMn^+ - PriiM+ne - Im+Tra Vj,j^j+npnj^j+, 

(A2) 
with um the metal number density, /3 — 3 x 10 ^"^ 
(100 K/T)^ " m^/s the dissociative recombination rate coef- 
ficient for molecular ions, /3^ = 3 x 10"^* (100 K/T')^/^ m^/s 
the radiative recombination rate coefficient for metal atoms 
and Pt = 3x 10~^^m^/s the rate coefficient of charge transfer 
from molecular ions to metal atoms. 

Equation ([3| is retrieved by summing these equations, 
recalling that rii = n^+ + 71^+ , if one puts: 



1 



•n-M+ 



^ 



+ ^m + 



V"^m + 



h 



Im+ ^A/+ I -'m+ ^m + 



V"T'M+ rii ym„+ Ui 



Pcff ~ Pr H P 



At equilibrium, provided that 

Prnlf-nt, 4 



XePrPsO. 



hi+ '"TM+ Tra'^npU'lf 3 ^1^+ Vtm+ "^^2 
that is, if we plug in equation Kt for Xe, 



«1, 



Psa\^ ^memM+ (Pr 



C-f) 



UM+mlj^ P 



«1, 



(A3) 



(A4) 



(A5) 



(A6) 



(A7) 



which, given that Pr <^ P, is essentially guaranteed by in- 
equality (Im, we draw from equation (A2|: 



nm+ ^ 3 Im+ yTM+ "^H"2 e 

n-A/+ 4 



psaPtXM 



10"^/ 



"+ (,100K 



1/2 



10-2 



10"^ kg/ m^ 
psa 



Xm 



APPENDIX B: LINEAR STABILITY ANALYSIS 

We consider the behavior of the equilibrium solution of our 
reduced model system ( |14[ ) and \1&\ to linear perturbations 
with a time dependence of the form exp (sf). The linearized 



form of equations ( 14 1 and (161 is 



'ps-=l,\p 



^'^+«l^^'-^(^)^^)) (^1) 



s5y^ = 2-,'{e)y^5e + 2 (7(e) - 2Ay^) 5y\ (B2) 

Using background equation (fTsl, equation ( |B2[ ) yields: 

i{e)y^5e 



Sy 



7(e) + S/2' 



(B3) 



Then, using equation ( |17[ ), and provided 27-l-s 7^ 0, equation 
(Bl I may be rewritten as: 



spSe = 



dz 



'¥s<'*: 



(B4) 



witiJEl 

h = 



7I7 + 



-jarctan 



i|7+f|iflms = 
2^±^^ifIms/0 ^^^^ 



2lo^Pv — '-ns7 
where Re and Im denote the real and imaginary parts, re- 



spectively. If we mutliply equation (B4| by [h&e)* and inte- 
grate over 2, we obtain: 



/- 



+ 00 Qy^ 

h 



Ji ihSe) \^pAz 



/+^ft*|5ePpd^ 



(B6) 



We now want to show that this equation implies that the 
real part of s (i.e. the growth rate of the perturbation) is 
negative, i.e. that the system is linearly stable. We distin- 
guish two cases: 

Case 1: Ims = 0. In this case h is real and positive. It 



is then clear from equation ( B6 1 that s must be real and 
negativsrl 

Case 2: Ims 7^ 0. We proceed ab absurdo, by supposing 
Res > 0. Then, 27 -I- Res is always strictly positive, from 



(A8) 



where xm = UM/nn- It is normalized to a reasonable 
value considering the solar abundances for e.g. Na, Mg, Si, 
K, Fe (log( M/i/) = -5.7 0, -4.45, -4.46, -6.89, -4.53, re- 
spectively ; Lodders 2003 1 and the expected 1-2 order-of- 



magnitude depletion due to condensation (see Section 2.1). 
Then, to a good approximation (if that ratio is larger than 
10~^), PcH ~ {nj-,-^+/nM+)P, justifying the normalization 



^ We use the fact that 

d / , , t- Reto \ 1 

— I ln|t — w] + jarctan I = , 

dt \ Imw J t — w 

for t a real variable and w a com plex constant. 
^ Recall that in writing equation | |B6| , we had assumed that 27+s 
was nonzero everywhere; if that is not the case, it is clear that s 
would be a negative real number in that case too. 
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which it follows, from the definition of h (equation (B5l) 
that: 



— f < arg/i < if Ims > 
< arg/i < f if Ims < 



,^..^, _..__. (B7) 



where argfe is the real number (between —n and tt) satisfying 
h — |/ijexp(iarg/i). In the complex plane, for a given eigen- 
vector (and thus a fixed Ims), h is thus confined to a deter- 
mined quadrant. Thus, the integrand in each of the integrals 



of ( B6 1 is always in the symmetric quadrant with respect to 
the real axis (because arg(l//i) = arg(/i*) — — arg/i). Same 
holds for the integrals themselves. It follows that the phase 
of the ratio of these integrals, which equals the difference of 
the phases of these two integrals, cannot be larger than 7r/2 
in absolute value. In other words, the ratio of the integrals 
has a positive real part and hence Res < 0. This contradicts 
the hypothesis and s must thus have a negative real part. 



